MODULE md1_LIP2_hermite

   USE md1_global_types
   USE md1_global_consts
   IMPLICIT NONE
   PRIVATE
   ! Public procedures
   PUBLIC :: LIP2_M1,     &
             LIP2_M2,     &
             LIP2_M3

CONTAINS

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M1(t,p,Q,lambda)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M1

      SELECT CASE (t)
         CASE (0)
            LIP2_M1 = LIP2_M10(p,Q,lambda)
         CASE (1)
            LIP2_M1 = LIP2_M11(p,Q,lambda)
         CASE (2)
            LIP2_M1 = LIP2_M12(p,Q,lambda)
!         CASE (3)
!            LIP2_M1 = LIP2_M13(p,Q,lambda)
!         CASE (4)
!            LIP2_M1 = LIP2_M14(p,Q,lambda)
         CASE DEFAULT
            CALL fmm_quit('NYI for LIP2/AO!')
      END SELECT

   END FUNCTION LIP2_M1

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M2(t,p,Q,lambda)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M2

      SELECT CASE (t)
         CASE (0)
            LIP2_M2 = LIP2_M20(p,Q,lambda)
         CASE (1)
            LIP2_M2 = LIP2_M21(p,Q,lambda)
         CASE (2)
            LIP2_M2 = LIP2_M22(p,Q,lambda)
!         CASE (3)
!            LIP2_M2 = LIP2_M23(p,Q,lambda)
!         CASE (4)
!            LIP2_M2 = LIP2_M24(p,Q,lambda)
         CASE DEFAULT
            CALL fmm_quit('NYI for LIP2/AO!')
      END SELECT

   END FUNCTION LIP2_M2

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M3(t,p,Q,lambda)

      IMPLICIT NONE
      INTEGER(INTK), INTENT(IN) :: t
      REAL(REALK),   INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M3

      SELECT CASE (t)
         CASE (0)
            LIP2_M3 = LIP2_M30(p,Q,lambda)
         CASE (1)
            LIP2_M3 = LIP2_M31(p,Q,lambda)
         CASE (2)
            LIP2_M3 = LIP2_M32(p,Q,lambda)
!         CASE (3)
!            LIP2_M3 = LIP2_M33(p,Q,lambda)
!         CASE (4)
!            LIP2_M3 = LIP2_M34(p,Q,lambda)
         CASE DEFAULT
            CALL fmm_quit('NYI for LIP2/AO!')
      END SELECT

   END FUNCTION LIP2_M3

!-------------------------------------------------------------------------------
!
! AngMom = 0
!
!-------------------------------------------------------------------------------

   FUNCTION LIP2_M10(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M10

      LIP2_M10 =   &
     dble(p ** (-0.3D1 / 0.2D1)) * (0.2D1 * exp(dble(4 * p * lambda ** 2))  &
     * Q * sqrt(dble(p)) + 0.2D1 * exp(dble(p) * (dble(4 * lambda ** 2)     &
     + Q ** 2)) * Q ** 2 * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q)    &
     * dble(p) + exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))            &
     * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q) - 0.6D1                &
     * exp(dble(4 * p * lambda ** 2)) * dble(lambda) * sqrt(dble(p))        &
     - 0.6D1 * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))              &
     * dble(lambda) * Q * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q)     &
     * dble(p) + 0.4D1 * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))    &
     * dble(lambda ** 2) * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q)    &
     * dble(p) + 0.2D1 * exp(0.4D1 * dble(p) * Q * dble(lambda))            &
     * dble(lambda) * sqrt(dble(p)) - 0.2D1 * Q * exp(0.4D1 * dble(p) * Q   &
     * dble(lambda)) * sqrt(dble(p)) + 0.2D1 * dble(p) * Q ** 2             &
     * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda) * sqrt(dble(p))      &
     - sqrt(dble(p)) * Q) * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2)) &
     + sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda) * sqrt(dble(p))      &
     - sqrt(dble(p)) * Q) * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2)) &
     - 0.6D1 * dble(lambda) * dble(p) * Q * sqrt(0.3141592654D1)            &
     * erf(0.2D1 * dble(lambda) * sqrt(dble(p)) - sqrt(dble(p)) * Q)        &
     * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2)) + 0.4D1              &
     * dble(lambda ** 2) * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)  &
     * sqrt(dble(p)) - sqrt(dble(p)) * Q) * dble(p) * exp(dble(p)           &
     * (dble(4 * lambda ** 2) + Q ** 2))) * exp(-dble(p)                    &
     * (dble(4 * lambda ** 2) + Q ** 2)) / dble(lambda ** 2) / 0.8D1

   END FUNCTION LIP2_M10

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M20(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M20

      LIP2_M20 =   &
        dble(p ** (-0.3D1 / 0.2D1)) * (-0.2D1 * exp(dble(4*p*lambda**2))    &
      * Q * sqrt(dble(p)) - 0.2D1 * exp(dble(p) * (dble(4 * lambda ** 2)    &
      + Q ** 2)) * Q ** 2 * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q)   &
      * dble(p) - exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))           &
      * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q) + 0.4D1 * exp(dble(4  &
      * p * lambda ** 2)) * dble(lambda) * sqrt(dble(p)) + 0.4D1            &
      * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2)) * dble(lambda)      &
      * Q * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q) * dble(p) + 0.2D1 &
      * Q * exp(0.4D1 * dble(p) * Q * dble(lambda)) * sqrt(dble(p)) - 0.2D1 &
      * dble(p) * Q ** 2 * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)  &
      * sqrt(dble(p)) - sqrt(dble(p))*Q) * exp(dble(p) * (dble(4*lambda**2) &
      + Q ** 2)) - sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)          &
      * sqrt(dble(p)) - sqrt(dble(p)) *Q) *exp(dble(p) * (dble(4*lambda**2) & 
      + Q ** 2)) + 0.4D1 * dble(lambda) * dble(p) * Q *sqrt(0.3141592654D1) &
      * erf(0.2D1 * dble(lambda) * sqrt(dble(p)) - sqrt(dble(p)) * Q)       &
      * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))) * exp(-dble(p)     &
      * (dble(4 * lambda ** 2) + Q ** 2)) / dble(lambda ** 2) / 0.4D1

   END FUNCTION LIP2_M20

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M30(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M30

      LIP2_M30 =   &
     -dble(p ** (-0.3D1 / 0.2D1)) * (-0.2D1 * exp(dble(4 * p * lambda ** 2)) &
     * Q * sqrt(dble(p)) - 0.2D1 * exp(dble(p) * (dble(4 * lambda ** 2)      &
     + Q ** 2)) * Q ** 2 * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q)     &
     * dble(p) - exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))             &
     * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q) + 0.2D1                 &
     * exp(dble(4 * p * lambda ** 2)) * dble(lambda) * sqrt(dble(p)) + 0.2D1 &
     * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2)) * dble(lambda) * Q    &
     * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * Q) * dble(p) + 0.2D1       &
     * exp(0.4D1 * dble(p) * Q * dble(lambda)) * dble(lambda) *sqrt(dble(p)) &
     + 0.2D1 * Q * exp(0.4D1 * dble(p) * Q * dble(lambda)) * sqrt(dble(p))   &
     - 0.2D1 * dble(p) * Q ** 2 * sqrt(0.3141592654D1) * erf(0.2D1           &
     * dble(lambda) * sqrt(dble(p)) - sqrt(dble(p)) * Q) * exp(dble(p)       &
     * (dble(4 * lambda ** 2) + Q ** 2)) - sqrt(0.3141592654D1) * erf(0.2D1  &
     * dble(lambda) * sqrt(dble(p)) - sqrt(dble(p)) * Q) * exp(dble(p)       &
     * (dble(4 * lambda ** 2) + Q ** 2)) + 0.2D1 * dble(lambda) * dble(p)    &
     * Q * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda) * sqrt(dble(p))   &
     - sqrt(dble(p)) * Q) * exp(dble(p) * (dble(4 * lambda ** 2) + Q ** 2))) & 
     * exp(-dble(p) * (dble(4*lambda**2) + Q**2)) / dble(lambda ** 2) / 0.8D1

   END FUNCTION LIP2_M30

!-------------------------------------------------------------------------------
!
! AngMom = 1
!
!-------------------------------------------------------------------------------

   FUNCTION LIP2_M11(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M11

      LIP2_M11 = -exp(-dble(p * (4 * lambda ** 2 + Q ** 2))) * (0.3D1    &
        * exp(dble(p * (4 * lambda ** 2 + Q ** 2))) * dble(lambda)    &
        * sqrt(0.3141592654D1) * erf(dble(Q) * sqrt(dble(p))) * dble(p)    &
        + 0.2D1 * exp(dble(4 * p * Q * lambda)) * sqrt(dble(p)) - 0.2D1    &
        * exp(dble(4 * p * lambda ** 2)) * sqrt(dble(p)) - 0.2D1    &
        * exp(dble(p * (4 * lambda ** 2 + Q ** 2))) * dble(Q)    &
        * sqrt(0.3141592654D1) * erf(dble(Q) * sqrt(dble(p))) * dble(p)    &
        - 0.2D1 * dble(p) * dble(Q) * sqrt(0.3141592654D1) * erf(0.2D1    &
        * sqrt(dble(p)) * dble(lambda) - dble(Q) * sqrt(dble(p)))    &
        * exp(dble(p * (4 * lambda ** 2 + Q ** 2))) - 0.4D1 * exp(dble(4    &
        * p * lambda ** 2)) * dble(lambda ** 2) * dble(p ** (0.3D1    &
        / 0.2D1)) + 0.3D1 * dble(lambda) * dble(p) * sqrt(0.3141592654D1)    &
        * erf(0.2D1 * sqrt(dble(p)) * dble(lambda) - dble(Q)    &
        * sqrt(dble(p))) * exp(dble(p * (4 * lambda ** 2 + Q ** 2))))    &
        * dble(p ** (-0.3D1 / 0.2D1)) / dble(lambda ** 2) / 0.4D1

   END FUNCTION LIP2_M11

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M21(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M21

      LIP2_M21 = exp(-dble(p * (4 * lambda ** 2 + Q ** 2)))    &
        * (-exp(dble(p * (4 * lambda ** 2 + Q ** 2))) * dble(Q)    &
        * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * dble(Q)) * dble(p)    &
        - exp(dble(4 * p * lambda ** 2)) * sqrt(dble(p))    &
        + exp(dble(4 * p * Q * lambda)) * sqrt(dble(p))    &
        + exp(dble(p * (4 * lambda ** 2 + Q ** 2))) * dble(lambda)    &
        * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * dble(Q)) * dble(p)    &
        - dble(p) * dble(Q) * sqrt(0.3141592654D1) * erf(0.2D1    &
        * dble(lambda) * sqrt(dble(p)) - sqrt(dble(p)) * dble(Q))    &
        * exp(dble(p * (4 * lambda ** 2 + Q ** 2))) + dble(lambda)    &
        * dble(p) * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)    &
        * sqrt(dble(p)) - sqrt(dble(p)) * dble(Q)) * exp(dble(p    &
        * (4 * lambda ** 2 + Q ** 2)))) * dble(p ** (-0.3D1 / 0.2D1))    &
        / dble(lambda ** 2)

   END FUNCTION LIP2_M21

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M31(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M31

      LIP2_M31 = -exp(-dble(p * (4 * lambda ** 2 + Q ** 2))) * (-0.2D1   &
        * exp(dble(4 * p * lambda ** 2)) * sqrt(dble(p)) + 0.4D1   &
        * dble(p ** (0.3D1 / 0.2D1)) * dble(lambda ** 2)   &
        * exp(dble(4 * p * Q * lambda)) + 0.2D1 * exp(dble(4 * p   &
        * Q * lambda)) * sqrt(dble(p)) - 0.2D1 * exp(dble(p * (4   &
        * lambda ** 2 + Q ** 2))) * dble(Q) * sqrt(0.3141592654D1)   &
        * erf(sqrt(dble(p)) * dble(Q)) * dble(p) + exp(dble(p * (4   &
        * lambda ** 2 + Q ** 2))) * dble(lambda) * sqrt(0.3141592654D1)   &
        * erf(sqrt(dble(p)) * dble(Q)) * dble(p) - 0.2D1 * dble(p)   &
        * dble(Q) * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)   &
        * sqrt(dble(p)) - sqrt(dble(p)) * dble(Q)) * exp(dble(p   &
        * (4 * lambda ** 2 + Q ** 2))) + dble(lambda) * dble(p)   &
        * sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda) * sqrt(dble(p))   &
        - sqrt(dble(p)) * dble(Q)) * exp(dble(p * (4 * lambda ** 2   &
        + Q ** 2)))) * dble(p ** (-0.3D1 / 0.2D1))   &
        / dble(lambda ** 2) / 0.4D1

   END FUNCTION LIP2_M31

!-------------------------------------------------------------------------------
!
! AngMom = 2
!
!-------------------------------------------------------------------------------

   FUNCTION LIP2_M12(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M12

      LIP2_M12 = -exp(-dble(p * (4 * lambda ** 2 + Q ** 2))) * (0.4D1   &
        * dble(Q) * dble(p ** (0.3D1 / 0.2D1)) * exp(dble(4 * p   &
        * lambda ** 2)) * dble(lambda ** 2) + 0.3D1 * exp(dble(4 * p   &
        * lambda ** 2)) * dble(lambda) * sqrt(dble(p)) + exp(dble(4 * p   &
        * Q * lambda)) * dble(lambda) * sqrt(dble(p)) - exp(dble(p * (4   &
        * lambda ** 2 + Q ** 2))) * sqrt(0.3141592654D1) * erf(dble(Q)   &
        * sqrt(dble(p))) - sqrt(0.3141592654D1) * erf(0.2D1 * sqrt(dble(p))   &
        * dble(lambda) - dble(Q) * sqrt(dble(p))) * exp(dble(p * (4   &
        * lambda ** 2 + Q ** 2)))) * dble(p ** (-0.1D1 / 0.2D1))   &
        / dble(lambda ** 2) / 0.2D1

   END FUNCTION LIP2_M12

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M22(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M22

      LIP2_M22 = exp(-dble(p * (4 * lambda ** 2 + Q ** 2)))    &
        * (-exp(dble(p * (4 * lambda ** 2 + Q ** 2)))    &
        * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * dble(Q)) + 0.2D1    &
        * sqrt(dble(p)) * dble(lambda) * exp(dble(4 * p * Q * lambda))    &
        + 0.2D1 * exp(dble(4 * p * lambda ** 2)) * dble(lambda)    &
        * sqrt(dble(p)) - sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda)    &
        * sqrt(dble(p)) - sqrt(dble(p)) * dble(Q)) * exp(dble(p * (4    &
        * lambda ** 2 + Q ** 2)))) * dble(p ** (-0.1D1 / 0.2D1))    &
        / dble(lambda ** 2)

   END FUNCTION LIP2_M22

!-------------------------------------------------------------------------------

   FUNCTION LIP2_M32(p,Q,lambda)

      IMPLICIT NONE
      REAL(REALK), INTENT(IN) :: p, Q, lambda
      REAL(REALK) :: LIP2_M32

      LIP2_M32  = -exp(-dble(p * (4 * lambda ** 2 + Q ** 2))) * (-0.4D1   &
        * dble(Q) * dble(p ** (0.3D1 / 0.2D1)) * dble(lambda ** 2)   &
        * exp(dble(4 * p * Q * lambda)) + 0.8D1 * dble(p ** (0.3D1   &
        / 0.2D1)) * dble(lambda ** 3) * exp(dble(4 * p * Q * lambda))   &
        + 0.3D1 * exp(dble(4 * p * Q * lambda)) * dble(lambda)   &
        * sqrt(dble(p)) - exp(dble(p * (4 * lambda ** 2 + Q ** 2)))   &
        * sqrt(0.3141592654D1) * erf(sqrt(dble(p)) * dble(Q))   &
        + exp(dble(4 * p * lambda ** 2)) * dble(lambda) * sqrt(dble(p))   &
        - sqrt(0.3141592654D1) * erf(0.2D1 * dble(lambda) * sqrt(dble(p))   &
        - sqrt(dble(p)) * dble(Q)) * exp(dble(p * (4 * lambda ** 2   &
        + Q ** 2)))) * dble(p ** (-0.1D1 / 0.2D1)) / dble(lambda ** 2)  &
        / 0.2D1

   END FUNCTION LIP2_M32

!-------------------------------------------------------------------------------

END MODULE md1_LIP2_hermite

